function [thetaFMRate] = rvdiffeq(thetaFM, omegaFMM)
%RVDIFFEQ 等效旋转矢量微分方程（Bortz方程）
%   此处显示详细说明

theta = norm(thetaFM);
thetaCrossomega = cross(thetaFM, omegaFMM);
if isnumeric(thetaFM) && isnumeric(omegaFMM) && theta<0.05
    % when theta is smaller than 3 degrees, use taylor expansion to avoid singularity
    coef = theta^4/30240 + theta^2/720 + 1/12;
else
    coef = (1-theta*sin(theta)/2/(1-cos(theta)))/theta^2;
end
thetaFMRate = omegaFMM + thetaCrossomega/2 + coef*cross(thetaFM, thetaCrossomega);
end